function [L,dL] = liklihoodFun(q,path)



L = 0;
% dL = sym('A%d%d',[8,8]);
% dL(:,:) = 0;
dL = zeros(0,0);
DEFAULT = 8;
NUM_RATE = 7;
MIN_TIME_UNIT = 1/12; % 1 month in year measure


for i = 1: length(path)
    temp = path(i).totalsRec;
    % check if the rating is same
    if temp(1,2)~= temp(1,3);
        % for cases with rating changes
        if isempty(find(temp == DEFAULT))
            % censoring situation
            length_temp = size(temp,1);
            for k = 1: length_temp
                if k ~= length_temp
                    p = expm(q*temp(k,1));
                    L = L + log(p(temp(k,2),temp(k,3)));
                    dL= dL+ derivativeL(NUM_RATE,temp(k,1),temp(k,2), temp(k,3), q, 'normal');
                    
                    
                else
                    % deal with censoring case
                    Pij = 0;
                    for l = 1:NUM_RATE
                        p = expm(q*temp(k,1));
                        Pij = Pij + p(temp(k,3),l);
                    end
                    
                    L = L + log(Pij);
                    dL= dL+ derivativeL(NUM_RATE,temp(k,1),temp(k,2), temp(k,3), q, 'cencor');
                    
                end
            end
            
        else % default situation
            [ROW,COL] = find(temp ==DEFAULT);
            ROW_MIN = min(ROW);
            
            for k = 1: ROW_MIN
                if k ~= ROW_MIN
                    p = expm(q*temp(k,1));
                    L = L+ log( p(temp(k,2),temp(k,3)) );
                    dL= dL+ derivativeL(NUM_RATE,temp(k,1),temp(k,2), temp(k,3), q, 'normal');
                    
                else
                    % deal with default case
                    Pij = 0;
                    t1 = (temp(k,1)-MIN_TIME_UNIT);
                    t2 =  MIN_TIME_UNIT;
                    p1 = expm(q*t1);
                    p2 = expm(q*t2);
                    for l = 1:NUM_RATE
                        
                        a = p1(temp(k,2),l);
                        b = p2(l,temp(k,3));
                        Pij = Pij + a*b;
                    end
                    
                    L = L + log(Pij);
                    dL= dL+ derivativeL(NUM_RATE,temp(k,1),temp(k,2), temp(k,3), q, 'default');
                
                end
                
            end
            
        end
        
    end
end







end


%                     for para_i = 1: NUM_RATE
%                         for para_j = 1:NUM_RATE+1
%                             if para_i ~= para_j
%                                 T = temp(k,1);
%                                 dQ = derivativeQ(para_i, para_j);
%                                 dP = derivativeP(T,dQ,q);
%                                 p = expm(q*T);
%                                 dL(para_i,para_j) = dL(para_i,para_j) + 1/p(temp(k,2),temp(k,3))* dP(temp(k,2),temp(k,3));
%                             end
%                         end
%                     end



% function p = prob(q,t)
% [A,D] = eig(q);
% p = A*diag(exp(diag(t*D)))/A;% D* Inv(A) ==> D/A
% end